#include "main_control.h"
#include "nvtx3/nvToolsExt.h"

void Main_Control::Step_I() {
  Time_Count TC;
  Time_Count TC_Tot;
  TC.start();
  TC_Tot.start();
  this->print_General();
  // Get Input
  cout << " +++  Step_I +++   " << endl;
  get_info.init(Filename_Input);
  get_info.build();
  string f_dst = "serial_dat/Output/DF.Input";
  File_M.copy(Filename_Input, f_dst);
  // Init coeff
  nvtxRangePush("init_coeff");
  this->init_coeff();
  nvtxRangePop();
  nvtxRangePush("build_basis");
  this->build_basis();
  nvtxRangePop();
  cout << " -- Finished State & Coeff successfully  -- " << endl;
  //
  int Trans_Test_Flag = -1;
  Trans_Test_Flag = get_info.trans_test.Trans_Test_Flag;
  if (Trans_Test_Flag == 1) {
    TC.start();
    cout << "T-coefficients Start : " << endl;
    Hobra hobra_t(bin_t, cg_x, ninJ_x);
    nvtxRangePush("Trans_coeff");
    TransCoeff trans_t(sixJ_x, ninJ_x, hobra_t, jaco_nl, single_t, triple_t);
    //trans_t.build();
    trans_t.build_OA();
    nvtxRangePop();
    string trans_file = "serial_dat/Output/trans_vec.dat";
    string trans_file_ser = "serial_dat/Output/trans_serial.dat";
    // if (Print_Mid_Flag > 0) {
    // ser_t.save(trans_t, trans_file_ser);
    trans_t.print(trans_file);
    //}
    cout << "T-coefficients Time cost : " << TC.sec() << "s" << endl;
    Check_File CF;
    int flag_diff = -1;
    string trans_file_std = get_info.trans_test.trans_file_std;
    flag_diff = CF.diff_two_file(trans_file, trans_file_std);
    if (flag_diff == 1) {
      cout << " ==  Correct Trans ! == " << endl;
    }
    cout << " Step I total Time cost : " << TC_Tot.min() << " min  or "
         << TC_Tot.hour() << " h on " << MainCom.size() << "nodes" << endl;

    exit(0);
  }
  // --- construct Wave ---
  vector<vector<int>> E_Jac_conf;
  E_Jac_conf = triple_t.E_Jac_conf;
  const int lamb_cut = get_info.wave_info.Lambda_Cut;
  double hbar_omg = get_info.wave_info.hbar_omega;
  int reg_pow = get_info.wave_info.regulator_power;
  int k_mesh = get_info.wave_info.k_mesh;
  int k_scale = get_info.wave_info.k_scale;
  int p_mesh = get_info.wave_info.p_mesh;
  int Emax_F = get_info.wave_info.EF_max;
  int Emax_I = get_info.wave_info.EI_max;
  int n_max = get_info.wave_info.n_max;
  int l_max = get_info.wave_info.l_max;
  int lbarS = get_info.wave_info.lbarS;
  int lbarS_G = get_info.wave_info.lbarS_G;
  int lbarS_P = get_info.wave_info.lbarS_P;

  int EJC_num = E_Jac_conf.size();
  int EJC_max = 0;
  for (int i = 0; i < EJC_num; i++) {
    if (E_Jac_conf[i][0] > EJC_max)
      EJC_max = E_Jac_conf[i][0];
    if (E_Jac_conf[i][1] > EJC_max)
      EJC_max = E_Jac_conf[i][1];
  }
  string wave_serial_file =
      "../Wave_mpi_E/serial_dat/Wave_serial" + std::to_string(lamb_cut) +
      "_ho" + std::to_string(hbar_omg) + "_reg" + std::to_string(reg_pow) +
      "_k" + std::to_string(k_mesh) + "_ks" + std::to_string(k_scale) + "_p" +
      std::to_string(p_mesh) + "_EF" + std::to_string(Emax_F) + "_EI" +
      std::to_string(Emax_I) + "_EJC" + std::to_string(EJC_num) + "-" +
      std::to_string(EJC_max) + "_n" + std::to_string(n_max) + "_l" +
      std::to_string(l_max) + "_g" + std::to_string(lbarS) + "_G" +
      std::to_string(lbarS_G) + "_P" + std::to_string(lbarS_P) + "_sep.dat";
  cout << wave_serial_file << endl;
  Wave_Fun wave(lamb_cut, hbar_omg, reg_pow, k_mesh, k_scale, p_mesh, Emax_F,
                Emax_I);
  cout << "Test press :  Yes " << endl;

  // --- Press Test ---
  cout << "TP all_size = " << triple_t.all_size_T << endl;
  cout << "Jacobi state all_size = " << jaco_b.all_size;
  cout << "\t Jacobi state Emax = " << jaco_b.Emax << endl;
  cout << "\t E_tot : " << endl;
  for (int k = 0; k < triple_t.TPconf_E_tot.size(); k++) {
    cout << "File No : " << k << endl;
    for (int i = 0; i < triple_t.TPconf_E_tot[k].size(); i++) {
      cout << "\t" << triple_t.TPconf_E_tot[k][i].E_F;
      cout << "\t" << triple_t.TPconf_E_tot[k][i].E_I << endl;
    }
  }
  cout << "----- E_cut_conf ------ " << endl;
  vector<vector<int>> E_cut_conf;
  triple_t.find_E_Jac_conf(E_cut_conf);
  cout << E_cut_conf.size() << "\t E_conf Num" << endl;
  int EF_max = -1;
  int EI_max = -1;
  for (int i = 0; i < E_cut_conf.size(); i++) {
    cout << E_cut_conf[i][0] << "\t " << E_cut_conf[i][1] << "\t E_F \t E_I"
         << endl;
    if (E_cut_conf[i][0] > EF_max)
      EF_max = E_cut_conf[i][0];
    if (E_cut_conf[i][1] > EI_max)
      EI_max = E_cut_conf[i][1];
  }
  cout << "EF_max = " << EF_max << "\t EI_max = " << EI_max << endl;
  get_info.wave_info.EF_max = EF_max;
  get_info.wave_info.EI_max = EI_max;
  cout << "Jacobi_NL state all_size = " << jaco_nl.all_size << endl;
  int size_trans = 0;
  /*for (int i = 0; i < trans_t.Tval_vec.size(); i++) {
    for (int j = 0; j < trans_t.Tval_vec[i].size(); j++) {
      for (int k = 0; k < trans_t.Tval_vec[i][j].size(); k++) {
        size_trans += trans_t.Tval_vec[i][j][k].size();
      }
    }
  } */
  cout << " Trans-coeff size = " << size_trans << "\t with Memory ~ ";
  cout << size_trans * 8.0 / 1024.0 / 1024.0 << " Mb" << endl;

  cout << "nmax = " << jaco_b.nmax << "\t lmax = " << jaco_b.lmax
       << "\t J12_2max = " << jaco_b.J12_2max
       << "\t j3_2max = " << jaco_b.j3_2max << endl;
  get_info.wave_info.n_max = jaco_b.nmax;
  get_info.wave_info.l_max = jaco_b.lmax;
  Tpe_term tpe_test(LECs, wave, cg_x, sixJ_x, ninJ_x, jaco_b);
  // cout << "tpe_test start to init" << endl;
  tpe_test.init(jaco_b.Jaco_l123_info, jaco_b.Jaco_E_info, jaco_b.J12_2max,
                jaco_b.j3_2max);
  tpe_test.get_Llbar_info();
  cout << "Ic1 size = " << tpe_test.Ic1_size << endl;
  cout << "Ic3 size = " << tpe_test.Ic3_size << endl;
  //<< "\t Ic3 lbarS = " << tpe_test.Ic3_lbar_max << endl;
  cout << "Ic4 size = " << tpe_test.Ic4_size
       << "\t Ic4 lbarS = " << tpe_test.Ic4_lbar_max << endl;
  int lbar_MAX = tpe_test.Ic4_lbar_max;
  cout << " lbar_MAX : " << lbar_MAX << endl;
  get_info.wave_info.lbarS = lbar_MAX;
  get_info.wave_info.lbarS_G = lbar_MAX;
  get_info.wave_info.lbarS_P = lbar_MAX;
  Info_Wave IW_t;
  IW_t.n_max = jaco_b.nmax;
  IW_t.l_max = jaco_b.lmax;
  IW_t.lbar123_max = tpe_test.lbar123_max;
  IW_t.E_Jac_conf = E_Jac_conf;
  IW_t.Llbar_Ic_Reduced = tpe_test.Llbar_Ic_Reduced;
  IW_t.Llbar_Ic_Red_vec = tpe_test.Llbar_Ic_Red_vec;
  IW_t.l_vec = tpe_test.l_vec;
  IW_t.lbar_vec = tpe_test.lbar_vec;
  IW_t.l_state.clear();
  int size_l = tpe_test.l_state.size();
  for (int l = 0; l < size_l; l++) {
    l_struct_IW l_t;
    l_t.l12p = tpe_test.l_state[l].l12p;
    l_t.l3p = tpe_test.l_state[l].l3p;
    l_t.l12 = tpe_test.l_state[l].l12;
    l_t.l3 = tpe_test.l_state[l].l3;
    IW_t.l_state.push_back(l_t);
  }
  int size_lbar = tpe_test.lbar_state.size();
  IW_t.lbar_state.clear();
  for (int i = 0; i < size_lbar; i++) {
    lbar_struct_IW lb_t;
    lb_t.lbar = tpe_test.lbar_state[i].lbar;
    lb_t.lbar12 = tpe_test.lbar_state[i].lbar12;
    lb_t.lbar3 = tpe_test.lbar_state[i].lbar3;
    IW_t.lbar_state.push_back(lb_t);
  }
  Serial_SL_IW SL_IW;
  string serial_file_iw = "serial_dat/Output/info_wave.dat";
  SL_IW.save(IW_t, serial_file_iw);
  string Filename_Input_press = "data/Input/DF_press_test.Input";
  get_info.print(Filename_Input_press);
  cout << "  Generate data/Input/DF_press_test.Input " << endl;
  cout << " Step I total Time cost : " << TC_Tot.min() << " min  or "
       << TC_Tot.hour() << " h on " << MainCom.size() << "nodes" << endl;
  cout << " --- Finished Press Test --- " << endl;
  std::ofstream step_status("Step_Status.dat");
  step_status << "1  // Finished Step_I " << endl;
  step_status.close();
  exit(0);
}

void Main_Control::Step_II() {
  Time_Count TC;
  TC.start();
  Time_Count TC_Tot;
  if (MainCom.rank() == 0) {

    TC_Tot.start();
    this->print_General();
    // Get Input
    cout << " +++  Step_II +++   " << endl;
  }
  this->load_coeff();
  this->load_basis();
  get_info.init(Filename_Input);
  if (MainCom.rank() == 0) {
    get_info.build();
  } else {
    get_info.build(-1);
  }
  // --- construct Wave ---

  int lamb_cut = get_info.wave_info.Lambda_Cut;
  double hbar_omg = get_info.wave_info.hbar_omega;
  int reg_pow = get_info.wave_info.regulator_power;
  int k_mesh = get_info.wave_info.k_mesh;
  int k_scale = get_info.wave_info.k_scale;
  int p_mesh = get_info.wave_info.p_mesh;
  int Emax_F = get_info.wave_info.EF_max;
  int Emax_I = get_info.wave_info.EI_max;

  Wave_Fun wave(lamb_cut, hbar_omg, reg_pow, k_mesh, k_scale, p_mesh, Emax_F,
                Emax_I);
  //  --- Finished construct Wave ---
  TC.start();
  if (MainCom.rank() == 0) {
    cout << "\t ---- ---- Start Part : Wave & Tpe ---- ---- " << endl;
    // cout << "Wave Fun Initial Time cost : " << TC.sec() << "s" << endl;
    // cout << " Wave lbar_glbarS_max = " << wave.lbar_glbarS_max << endl;
    // cout << " Wave lbar_PFglbarS_max = " << wave.lbar_PFglbarS_max << endl;
  }
  int mpi_size = MainCom.size();
  int mpi_rank = MainCom.rank();
  // int tpe_init_flag = 1;
  vector<int> mpi_run_flag_gather(MainCom.size());
  int mpi_run_flag = 1;
  all_gather(MainCom, mpi_run_flag, mpi_run_flag_gather);
  //
  Tpe_term tpe_t(LECs, wave, cg_x, sixJ_x, ninJ_x, jaco_b);
  cout << "nmax = " << jaco_b.nmax << "\t lmax = " << jaco_b.lmax
       << "\t J12_2max = " << jaco_b.J12_2max
       << "\t j3_2max = " << jaco_b.j3_2max << endl;
  vector<vector<int>> E_Jac_conf;
  E_Jac_conf = triple_t.E_Jac_conf;
  tpe_t.init_II(jaco_b.Jaco_l123_info, jaco_b.Jaco_E_info, E_Jac_conf,
                jaco_b.J12_2max, jaco_b.j3_2max);
  tpe_t.TPE_Build_MPI(E_Jac_conf);
  //
  MainCom.barrier();
  all_gather(MainCom, mpi_run_flag, mpi_run_flag_gather);
  if (MainCom.rank() == 0) {
    tpe_t.print("serial_dat/Output/tpe.dat");
    ser_t.save(tpe_t, "serial_dat/Output/Tpe_serial.dat");
  }
  all_gather(MainCom, mpi_run_flag, mpi_run_flag_gather);
  cout << "Finished Wave & Tpe and Exit successfully" << endl;
  all_gather(MainCom, mpi_run_flag, mpi_run_flag_gather);
  if (MainCom.rank() == 0) {
    cout << " Step II total Time cost : " << TC_Tot.min() << " min  or "
         << TC_Tot.hour() << " h on " << MainCom.size() << "nodes" << endl;
  }
  std::ofstream step_status("Step_Status.dat");
  step_status << "2  // Finished Step_II " << endl;
  step_status.close();
  exit(0);
}

void Main_Control::Step_III() {
  Time_Count TC;
  Time_Count TC_Tot;
  TC.start();
  TC_Tot.start();
  this->print_General();
  // Get Input
  cout << " +++  Step_III +++   " << endl;
  this->load_coeff();
  this->load_basis();
  get_info.init(Filename_Input);
  get_info.build();
  // --- construct Wave ---

  int lamb_cut = get_info.wave_info.Lambda_Cut;
  double hbar_omg = get_info.wave_info.hbar_omega;
  int reg_pow = get_info.wave_info.regulator_power;
  int k_mesh = get_info.wave_info.k_mesh;
  int k_scale = get_info.wave_info.k_scale;
  int p_mesh = get_info.wave_info.p_mesh;
  int Emax_F = get_info.wave_info.EF_max;
  int Emax_I = get_info.wave_info.EI_max;

  Wave_Fun wave(lamb_cut, hbar_omg, reg_pow, k_mesh, k_scale, p_mesh, Emax_F,
                Emax_I);
  wave.init_single(jaco_b.nmax, jaco_b.lmax);
  //  ---    ---   ---    ---
  //  --- construct TPE ---
  Tpe_term tpe_t(LECs, wave, cg_x, sixJ_x, ninJ_x, jaco_b);
  ser_t.load(tpe_t, "serial_dat/Output/Tpe_serial.dat");
  //  ---    ---   ---    ---
  Hobra hobra_t(bin_t, cg_x, ninJ_x);
  //  --- construct Trans ---
  TC.start();
  TransCoeff trans_t(sixJ_x, ninJ_x, hobra_t, jaco_nl, single_t, triple_t);
  trans_t.build();
  if (Print_Mid_Flag > 0) {
    ser_t.save(trans_t, "serial_dat/Output/trans_serial.dat");
    trans_t.print("serial_dat/Output/trans_vec.dat");
  }
  cout << "T-coefficients Time cost : " << TC.sec() << "s" << endl;
  //  --- Contact Term ---
  TC.start();
  Cont_term cont_x(LECs, wave, sixJ_x, jaco_b);
  cont_x.init(jaco_b.nmax); //( n_max )
  cont_x.build();
  if (Print_Mid_Flag > 0) {
    cont_x.print("serial_dat/Output/cont_term.dat");
    ser_t.save(cont_x, "serial_dat/Output/Cont_serial.dat");
  }
  cout << "Contact term Time cost : " << TC.sec() << "s" << endl;
  //  --- OPE Term ---
  TC.start();
  Ope_term ope_x(LECs, wave, cg_x, sixJ_x, ninJ_x, jaco_b);
  ope_x.init();
  ope_x.build();
  if (Print_Mid_Flag > 0) {
    ope_x.print("serial_dat/Output/ope_term.dat");
    ser_t.save(ope_x, "serial_dat/Output/Ope_serial.dat");
  }
  cout << "OPE term Time cost : " << TC.sec() << "s" << endl;
  //  --- Package ---
  Antisym antisym_t(sixJ_x, ninJ_x, hobra_t, jaco_b);
  ser_t.load(antisym_t, "serial_dat/Output/anti_serial.dat");
  //
  Package pack(hbar_omg, triple_t, jaco_nl, jaco_b, trans_t, antisym_t, cont_x,
               ope_x, tpe_t);
  TC.start();
  pack.build_JacAvec();
  cout << "Pack build_JacAvec Time cost : " << TC.min() << "min" << endl;
  TC.start();
  pack.build();
  cout << "Pack Time cost : " << TC.min() << "min" << endl;
  TC.start();
  pack.print("serial_dat/Output/package.dat");
  cout << "Pack print Time cost : " << TC.sec() << "s" << endl;
  TC.start();
  ser_t.save(pack, "serial_dat/Output/Package_serial.dat");
  cout << "Pack save Time cost : " << TC.sec() << "s" << endl;
  TC.start();

  int num_t = get_info.tp_al_info.size();
  vector<vector<vector<int>>> TP_allow_vec_gather;
  vector<int> Same_particle_Flag_gather;
  for (int i = 0; i < num_t; i++) {
    vector<vector<int>> TP_allow_vec(2);
    for (int j = 0; j < 2; j++) {
      TP_allow_vec[j].resize(6);
      for (int k = 0; k < 6; k++) {
        TP_allow_vec[j][k] = get_info.tp_al_info[i].tp_al_vec[k][j];
      }
    }
    TP_allow_vec_gather.push_back(TP_allow_vec);
    Same_particle_Flag_gather.push_back(
        get_info.tp_al_info[i].Same_particle_Flag);
  }
  // ------- Output print -------- //
  std::string Pack_ALL_File_name = "serial_dat/Output/3BME_DF.dat";
  std::ofstream Pack_ALL_File;
  Pack_ALL_File.open(Pack_ALL_File_name);
  this->print_General(Pack_ALL_File);
  LEC_CONST LECs;
  Pack_ALL_File << " ---- Parameters : ----" << endl;
  Pack_ALL_File << "\t Lambda_Cut = " << get_info.wave_info.Lambda_Cut
                << "  (MeV)" << endl;
  Pack_ALL_File << "\t c1 = " << LECs.LEC_C1 << "  (GeV^-1)" << endl;
  Pack_ALL_File << "\t c3 = " << LECs.LEC_C3 << "  (GeV^-1)" << endl;
  Pack_ALL_File << "\t c4 = " << LECs.LEC_C4 << "  (GeV^-1)" << endl;
  Pack_ALL_File << "\t cD = " << LECs.LEC_CD << endl;
  Pack_ALL_File << "\t cE = " << LECs.LEC_CE << endl;
  Pack_ALL_File << "\t gA = " << LECs.LEC_gA << endl;
  Pack_ALL_File << "\t Hbar-Omega = " << get_info.wave_info.hbar_omega << endl;
  Pack_ALL_File << "\t The Power of Regulator : "
                << get_info.wave_info.regulator_power << endl;

  Pack_ALL_File << " ---- 3BME range ----  " << endl;
  // int num_t = Same_particle_Flag_g.size();
  for (int i = 0; i < num_t; i++) {
    int Same_particle_Flag = Same_particle_Flag_gather[i];
    if (Same_particle_Flag > 0) {
      Pack_ALL_File
          << "!! The index of <a'b'c'|V|abc> has the index constrain !!"
          << endl;
      Pack_ALL_File << "    ";
      if (Same_particle_Flag == 1)
        Pack_ALL_File << "\t a' = a " << endl;
      if (Same_particle_Flag == 2)
        Pack_ALL_File << "\t b' = b " << endl;
      if (Same_particle_Flag == 3)
        Pack_ALL_File << "\t c' = c " << endl;
      if (Same_particle_Flag == 12)
        Pack_ALL_File << "\t a' = a & b' = b" << endl;
      if (Same_particle_Flag == 13)
        Pack_ALL_File << "\t a' = a & c' = c" << endl;
      if (Same_particle_Flag == 23)
        Pack_ALL_File << "\t b' = b & c' = c" << endl;
      if (Same_particle_Flag == 123)
        Pack_ALL_File << "\t a' = a & b' = b & c' = c" << endl;
      single_t.print(TP_allow_vec_gather[i], Pack_ALL_File);
    } else {
      Pack_ALL_File << "TP constrain type : " << get_info.tp_info.Type << endl;
      if (get_info.tp_info.Type == "E_max")
        Pack_ALL_File << "  !!!  E_max max(E_a + E_b + E_c) <= "
                      << get_info.tp_info.E_max << "   !!!" << endl;
      single_t.print(Pack_ALL_File);
    }
  }
  pack.print_TP_ALL_f(Pack_ALL_File, single_t);
  Pack_ALL_File.close();
  cout << "print Pack_ALL_File Time cost : " << TC.min() << "min" << endl;
  cout << "  --- Generate 3BME succeffully in " << Pack_ALL_File_name
       << "  ---- " << endl;
  cout << " Step III total Time cost : " << TC_Tot.min() << " min  or "
       << TC_Tot.hour() << " h on " << MainCom.size() << "nodes" << endl;
  cout << " ==== ALL FINISHED ==== " << endl;
  std::ofstream step_status("Step_Status.dat");
  step_status << "3  // Finished Step_III " << endl;
  step_status.close();
}

void Main_Control::build_basis() {
  Time_Count TC;
  TC.start();
  // Build Single-Paticle states
  if (get_info.sp_info.Type == "Read") {
    single_t.build_T(get_info.sp_info.sp_vec);
  } else if (get_info.sp_info.Type == "E_max") {
    single_t.build_T_E(get_info.sp_info.E_max);
  } else {
    cerr << " Wrong @ build SP " << endl;
    cerr << "\t get_info.sp_info.Type : " << get_info.sp_info.Type << endl;
    exit(0);
  }
  single_t.build_Tz();
  single_t.print_Tz();
  single_t.print("serial_dat/Output/SP.dat");
  single_t.print_Tz("serial_dat/Output/SP_Tz.dat");
  ser_t.save(single_t, "serial_dat/Output/SP_serial.dat");
  // Build Triple-Paticle states
  // Get TP info
  vector<vector<int>> E_Jac_conf_t;
  vector<int> abc_vec(6);
  int TP_allow_file_num;
  // int TP_E_max;
  vector<vector<vector<int>>> TP_allow_gather;
  vector<int> E_cut_gather;
  vector<int> Same_Particles_gather;
  vector<vector<int>> E_sign_gather;
  TP_allow_file_num = get_info.tp_al_info.size();
  for (int i = 0; i < TP_allow_file_num; i++) {
    vector<vector<int>> TP_allow_vec(2);
    for (int j = 0; j < 2; j++) {
      TP_allow_vec[j].resize(6);
      for (int k = 0; k < 6; k++) {
        TP_allow_vec[j][k] = get_info.tp_al_info[i].tp_al_vec[k][j];
      }
    }
    vector<int> E_sign(6);
    for (int k = 0; k < 6; k++) {
      E_sign[k] = get_info.tp_al_info[i].tp_al_vec[k][2];
    }
    int Same_Particles_Flag = get_info.tp_al_info[i].Same_particle_Flag;
    int E_cut = get_info.tp_al_info[i].E_cut;
    TP_allow_gather.push_back(TP_allow_vec);
    E_cut_gather.push_back(E_cut);
    E_sign_gather.push_back(E_sign);
    Same_Particles_gather.push_back(Same_Particles_Flag);
    if (get_info.tp_info.Type == "abc") {
      single_t.print(TP_allow_vec);
    } else {
      single_t.print();
    }
  }
  // ---- build TP  ----
  cout << "Start build_T" << endl;
  if (get_info.tp_info.Type == "abc") {
    int a_min = get_info.tp_info.abc_vec[0];
    int a_max = get_info.tp_info.abc_vec[1];
    int b_min = get_info.tp_info.abc_vec[2];
    int b_max = get_info.tp_info.abc_vec[3];
    int c_min = get_info.tp_info.abc_vec[4];
    int c_max = get_info.tp_info.abc_vec[5];
    // triple_t.build_T(single_t, a_min, a_max, b_min, b_max, c_min, c_max);
    triple_t.init_Constrain(TP_allow_gather);
    triple_t.build_T_Constrain(single_t, a_min, a_max, b_min, b_max, c_min,
                               c_max);
  } else if (get_info.tp_info.Type == "E_max") {
    int Emax_t = get_info.tp_info.E_max;
    // triple_t.build_T(single_t, Emax_t);
    triple_t.init_Constrain(TP_allow_gather);
    triple_t.build_T_Constrain(single_t, Emax_t);
  } else {
    cerr << " Wrong @ build TP " << endl;
    cerr << "\t get_info.tp_info.Type : " << get_info.tp_info.Type << endl;
    exit(0);
  }
  //cout << "Start build_T test point 1" << endl;

  // triple_t.Select_FI_T(single_t, TP_allow_gather);

  // triple_t.print("serial_dat/Output/TP_before_delete.dat");
  triple_t.Delet_Sort_T();
  // cout << "Finished Delet_Sort_T" << endl;
  triple_t.Delet_conf_T(single_t, E_sign_gather, E_cut_gather,
                        Same_Particles_gather);
  // cout << "Finished Delet_conf_T" << endl;
  triple_t.build_conf_T(single_t, E_sign_gather, E_cut_gather,
                        Same_Particles_gather);
  // cout << "Finished build_conf_T" << endl;
  triple_t.find_E_Jac_conf(E_Jac_conf_t);
  //cout << "Start build_T test point 2" << endl;
  std::ofstream E_conf_file;
  E_conf_file.open("serial_dat/Output/E_conf.Input");
  E_conf_file << E_Jac_conf_t.size() << "\t E_conf Num" << endl;
  for (int i = 0; i < E_Jac_conf_t.size(); i++) {
    E_conf_file << E_Jac_conf_t[i][0] << "\t " << E_Jac_conf_t[i][1]
                << "\t E_F \t E_I" << endl;
  }
  E_conf_file.close();
  //cout << "Start build_T test point 3" << endl;
  triple_t.print_TPconf_info("serial_dat/Output/TPconf.dat",
                             "serial_dat/Output/TPconf_E.dat",
                             "serial_dat/Output/TPconf_E_tot.dat");
  triple_t.sort();
  triple_t.print("serial_dat/Output/TP.dat");
  triple_t.print_channel("serial_dat/Output/TP_channel.dat");
  // ---- build TP_Tz ----
  /*
  cout << "Start build_Tz" << endl;
  if (get_info.tp_info.Type == "abc") {
    int a_min = 2 * get_info.tp_info.abc_vec[0];
    int a_max = 2 * get_info.tp_info.abc_vec[1] + 1;
    int b_min = 2 * get_info.tp_info.abc_vec[2];
    int b_max = 2 * get_info.tp_info.abc_vec[3] + 1;
    int c_min = 2 * get_info.tp_info.abc_vec[4];
    int c_max = 2 * get_info.tp_info.abc_vec[5] + 1;
    triple_t.build_Tz(single_t, a_min, a_max, b_min, b_max, c_min, c_max);
  } else if (get_info.tp_info.Type == "E_max") {
    int Emax_t = get_info.tp_info.E_max;
    triple_t.build_Tz(single_t, Emax_t);
  } else {
    cerr << " Wrong @ build TP_Tz " << endl;
    cerr << "\t get_info.tp_info.Type : " << get_info.tp_info.Type << endl;
    exit(0);
  }
  triple_t.Select_FI_Tz(TP_allow_gather);
  triple_t.Delet_Sort_Tz();
  triple_t.Delet_conf_Tz(single_t, E_sign_gather, E_cut_gather,
                         Same_Particles_gather);
  triple_t.build_conf_Tz(single_t, E_sign_gather, E_cut_gather,
                         Same_Particles_gather);
  triple_t.print_TPconf_info_Tz("serial_dat/Output/TPconf_Tz.dat",
                                "serial_dat/Output/TPconf_E_Tz.dat");
  triple_t.print_Tz("data/TP_Tz.dat");
  triple_t.print_channel_Tz("data/TP_channel_Tz.dat");
  */
  ser_t.save(triple_t, "serial_dat/Output/TP_serial.dat");
  cout << "Finished build_T" << endl;
  //
  // ---- Jacobi ----
  cout << "Start Jacobi " << endl;
  int Emax_JC = triple_t.Emax;
  int J2max_JC=triple_t.J2_max_TP + 2 * triple_t.Emax;
  jaco_b.build_JT(Emax_JC, J2max_JC);
  //cout << "Start Jacobi test point 1" << endl;
  jaco_b.print("serial_dat/Output/Jaco_unselect.dat");
  jaco_b.select_FI_T(triple_t.TPconf_E_tot);
  jaco_b.Delet_Sort_T(triple_t);
  cout << "Start Jacobi test point 2" << endl;
  jaco_b.sort();
  jaco_b.print("serial_dat/Output/Jaco_unselect2.dat");
  jaco_b.print_channel("serial_dat/Output/Jaco_channel_unselect2.dat");
  // -- Antisym --
  cout << " Start Antisym " << endl;
  Hobra hobra_t(bin_t, cg_x, ninJ_x);
  Antisym antisym_t(sixJ_x, ninJ_x, hobra_t, jaco_b);
  cout << " Start Antisym test point 1  " << endl;
  double Anti_test =
      (1 - 2 * antisym_t.mat_cal(jaco_b.state[1][0], jaco_b.state[1][0])) / 3;
  cout << " \t ===+++=== TEST === Anti_test = " << Anti_test << endl;
  antisym_t.cal();
  antisym_t.print("serial_dat/Output/anti_vec.dat");
  jaco_b.print_channel("serial_dat/Output/Jaco_channel.dat");
  jaco_b.select();
  jaco_b.print("serial_dat/Output/Jaco_unselect3.dat");
  cout << "Jacobi State Finished & Antisymmetry Time cost : " << TC.sec() << "s"
       << endl;
  // --  Jacobi_NL  --
  cout << "Start to build Jaco_NL" << endl;
  cout << "\t Triple Emax = " << triple_t.Emax << endl;
  cout << "\t Triple Emax_Final = " << triple_t.Emax_Final;
  cout << "\t Triple Emax_Init = " << triple_t.Emax_Init << endl;
  jaco_nl.build_NL(jaco_b, triple_t);
  cout << "Finished build_NL" << endl;
  jaco_nl.select_FI_T(triple_t);
  cout << "Finished select_FI_T" << endl;
  vector<int> Jaco_used_Flag;
  jaco_nl.Delete_unused_Jaco_channel(jaco_b, Jaco_used_Flag);
  jaco_b.Delete_unused_Jaco_channel(Jaco_used_Flag);
  jaco_nl.sort();
  jaco_b.build_conf_T(triple_t.TPconf_E_tot);
  jaco_b.print_Jacconf("serial_dat/Output/Jacconf.dat");
  jaco_nl.print("serial_dat/Output/Jaco_NL.dat");
  jaco_nl.print_channel("serial_dat/Output/Jaco_NL_channel.dat");
  jaco_b.print("serial_dat/Output/Jaco.dat");
  jaco_b.print_channel("serial_dat/Output/Jaco_channel.dat");
  cout << "Jaco l12max = " << jaco_b.l12max << "\t l3max = " << jaco_b.l3max
       << "\t lmax = " << jaco_b.lmax << endl;
  cout << "Jaco n12max = " << jaco_b.n12max << "\t n3max = " << jaco_b.n3max
       << "\t nmax = " << jaco_b.nmax << endl;
  cout << "Jaco Emax_F = " << jaco_b.Emax_F
       << "\t Jaco Emax_I = " << jaco_b.Emax_I << endl;
  //
  ser_t.save(antisym_t, "serial_dat/Output/anti_serial.dat");
  ser_t.save(jaco_b, "serial_dat/Output/jaco_serial.dat");
  ser_t.save(jaco_nl, "serial_dat/Output/jacoNL_serial.dat");
}

void Main_Control::load_basis() {
  ser_t.load(single_t, "serial_dat/Output/SP_serial.dat");
  if (MainCom.rank() == 0)
    cout << "Finished Load single_t" << endl;
  ser_t.load(triple_t, "serial_dat/Output/TP_serial.dat");
  if (MainCom.rank() == 0)
    cout << "Finished Load triple_t" << endl;
  ser_t.load(jaco_nl, "serial_dat/Output/jacoNL_serial.dat");
  // ser_t.load(antisym_t, "serial_dat/Output/anti_serial.dat");
  ser_t.load(jaco_b, "serial_dat/Output/jaco_serial.dat");
}

void Main_Control::init_coeff() {
  int J_max_cg = get_info.coeff_para.J_max_cg;
  int J_max_sixJ = get_info.coeff_para.J_max_sixJ;
  int n_max = get_info.coeff_para.n_max_bino;
  cout << " init_coeff : J_max_cg = " << J_max_cg
       << "  J_max_sixJ = " << J_max_sixJ << "  n_max = " << n_max << endl;
  cg_x.init(J_max_cg);
  sixJ_x.init(J_max_sixJ);
  ninJ_x.init(5, sixJ_x);
  bin_t.init(n_max);
  ser_t.save(cg_x, "serial_dat/Output/cg_serial.dat");
  ser_t.save(sixJ_x, "serial_dat/Output/sixJ_serial.dat");
  ser_t.save(bin_t, "serial_dat/Output/binom_serial.dat");
  ser_t.save(ninJ_x, "serial_dat/Output/ninJ_serial.dat");
  // hobra_t.init(bin_t, cg_x, ninJ_x);
}

void Main_Control::load_coeff() {
  ser_t.load(cg_x, "serial_dat/Output/cg_serial.dat");
  ser_t.load(sixJ_x, "serial_dat/Output/sixJ_serial.dat");
  ser_t.load(bin_t, "serial_dat/Output/binom_serial.dat");
  ser_t.load(ninJ_x, "serial_dat/Output/ninJ_serial.dat");
  // hobra_t.init(bin_t, cg_x, ninJ_x);
}

void Main_Control::print_General() {
  cout << "  ==========================================================="
       << endl;
  cout << "      __        __           ___  __   __   ___  __   ___    "
       << endl;
  cout << "     |  \\  /\\  |__) |__/    |__  |__) /  \\ |__  /__`  |      "
       << endl;
  cout << "     |__/ /--\\ |  \\ |  \\    |    |  \\ \\__/ |___ .__/  |      "
       << endl;
  cout << "                                                             "
       << endl;
  cout << "                         version 2.1                         "
       << endl;
  cout << "                     by Yuanzhuo 2020-June                   "
       << endl;
  cout << "    ------------------------------------------------------   "
       << endl;
  cout << endl;
  Time_Count TC;
  TC.print_time();
  cout << endl;
}

void Main_Control::print_General(std::ofstream &file) {
  file << "  ==========================================================="
       << endl;
  file << "      __        __           ___  __   __   ___  __   ___    "
       << endl;
  file << "     |  \\  /\\  |__) |__/    |__  |__) /  \\ |__  /__`  |     "
       << endl;
  file << "     |__/ /--\\ |  \\ |  \\    |    |  \\ \\__/ |___ .__/  |    "
       << endl;
  file << "                                                             "
       << endl;
  file << "    ------------------------------------------------------   "
       << endl;
  file << endl;
  Time_Count TC;
  TC.print_time(file);
  file << endl;
}

void Main_Control::print_LECs() {
  cout << "\t ---- LECs : ----" << endl;
  cout << "\t c1 = " << LECs.LEC_C1 << endl;
  cout << "\t c3 = " << LECs.LEC_C3 << endl;
  cout << "\t c4 = " << LECs.LEC_C4 << endl;
  cout << "\t cD = " << LECs.LEC_CD << endl;
  cout << "\t cE = " << LECs.LEC_CE << endl;
  cout << "\t ----------------" << endl;
}
